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Abstract. We here analyse numerical simulations of su- 
personic, hypersonic and magnetohydrodynamic turbu- 
lence that is free to decay. Our goals are to understand 
the dynamics of the decay and the characteristic proper- 
ties of the shock waves produced. This will be useful for 
interpretation of observations of both motions in molecu- 
lar clouds and sources of non-thermal radiation. 

We find that decaying hypersonic turbulence possesses 
an exponential tail of fast shocks and an exponential de- 
cay in time, i.e. the number of shocks is proportional to 
t exp(— ktv) for shock velocity jump v and mean initial 
wavenumber k. In contrast to the velocity gradients, the 
velocity Probability Distribution Function remains Gaus- 
sian with a more complex decay law. 

The energy is dissipated not by fast shocks but by a 
large number of low Mach number shocks. The power loss 
peaks near a low-speed turn-over in an exponential dis- 
tribution. An analytical extension of the mapping closure 
technique is able to predict the basic decay features. Our 
analytic description of the distribution of shock strengths 
should prove useful for direct modeling of observable emis- 
sion. We note that an exponential distribution of shocks 
such as we find will, in general, generate very low excita- 
tion shock signatures. 
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- ISM: clouds - ISM: kinematics and dynamics 



1. Introduction 

Many structures we observe in the Universe have been 
shaped by fluid turbulence. In astronomy, we often ob- 
serve high speed turbulence driven by supersonic ordered 
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motions such as jets, supernova shocks and stellar winds 
(e.g. Franco & Carraminana 1999). Hypersonic speeds, 
with Mach numbers above 10, are commonly encountered. 
Clearly, to understand the structure, we require a theory 
for supersonic turbulence. Here, we concentrate on decay- 
ing turbulence, such as could be expected in the wakes 
of bow shocks, in the lobes of radio galaxies or following 
explosive events. Two motivating questions are: how fast 
does supersonic turbulence decay when not continuously 
replenished and how can we distinguish decaying turbu- 
lence from other dynamical forms? The first question has 
been answered through recent numerical simulations de- 
scribed below. The answer to the second question is sought 
here. We look for a deep understanding of the dynamics 
and physics which control decaying supersonic turbulence. 
From this, and a following study of driven turbulence, we 
can derive the analytical characteristics and the observa- 
tional signatures pertaining to supersonic turbulence. We 
caution that we specify to uniform three dimensional tur- 
bulence with an isothermal equation of state, an initially 
uniform magnetic field and periodic boundary conditions. 



Numerical studies of decaying supersonic turbulence 
in three dimensions have revealed a power-law decay of 
the energy in time following a short low-loss period (Mac 
Low et al. 1998; Stone et al. 1998). Simulations of decay- 
ing subsonic and incompressible turbulence show similar 
temporal behaviour (e.g. Galtier et al. 1997), as discussed 
by Mac Low et al. (1999). In the numerical experiments, 
random Gaussian velocity fields were generated with small 
wavenumber disturbances. Magnetic fields were included 
of various strengths. Mac Low (1999) concluded that the 
decay is so rapid under all conditions that the motions we 
observe in molecular clouds must be continuously driven. 
In this work, we analyse the Mach 5 simulations from Mac 
Low et al. (1998) as well as a new Mach 50 simulation. 
The hypersonic run should best illustrate the mechanisms 
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behind the development and evolution of the shock field, 
possibly revealing asymptotic solutions. 

The major goal is to derive the spectrum of shocks 
(the Shock Probability Distribution Function) generated 
by turbulence. Shocks are often responsible for detailed 
bright features, such as filamentary and sheet structures, 
within which particles are highly excited. An example of 
a region which appears to contain a chaotic mixture of 
shocks, termed a 'Supersonic Turbulent Reactor' is the 
DR 21 molecular hydrogen outflow, driven by a collimated 
wind from a high mass young star (Smith et al. 1998). 
The shock spectrum is related to the molecular excita- 
tion, with weak shocks being responsible for rotational 
excitation and strong shocks for vibrational excitation. 

Previous studies of compressible turbulence have con- 
centrated on the density and velocity structure of the cold 
gas rather than the shocks. Three dimensional subsonic 
and transonic simulations (e.g. Porter et al. 1994; Fal- 
garone et al. 1994), two dimensional supersonic motions 
(Vazquez- Semadeni 1994) as well as three dimensional 
supersonic turbulence have been discussed (Vazquez- 
Semadeni et al. 1996, Padoan et al. 1998). One attempts 
to describe and fit the density and velocity structures ob- 
served in molecular clouds. This is often appropriate for 
the interpretation of clouds since, although the Mach num- 
ber is still high, the shock speeds are too low to produce 
bright features. The simulations analysed here are also be- 
ing interpreted by Mac Low & Ossenkopf (2000) in terms 
of density structure. 

Despite a diversity of theory, and an increase in ana- 
lytical knowledge, a succinct understanding of turbulence 
has not been attained (see Lesieur 1997). Therefore, we 
need not apologise for not fully interpreting the results 
for the supersonic case. We do not look for a Kolmogorov- 
inspired theory for two reasons. First, fully developed tur- 
bulence becomes increasingly non-Gaussian towards small 
scales. These intermittency effects dominate the statistics 
of velocity jumps in supersonic turbulence. Second: the 
strong compressibility implies that a wavenumber analy- 
sis is irrelevant since the energy spectrum of a shock or of 
a system of shocks is simply k~ 2 (e.g. Gotoh & Kraich- 
nan 1993). Note that the Kolmogorov-like spectra found 
by Porter et al. (1994) appeared at late times when the 
flow is clearly subsonic (and also note that even the initial 
RMS Mach number was only unity, which rather stretches 
the definition of supersonic turbulence). 

We attempt here to construct a physical model to 
describe the non-Gaussian Probability Density Functions 
(PDFs) for the shock waves. We adapt the mapping clo- 
sure analysis, as applied to Navier-Stokes (incompressible) 
and Burgers (one dimensional and pressure free) turbu- 
lence (Kraichnan 1990), to compressible turbulence. Phe- 
nomenological approaches, such as multifractal models or 
the log-normal model, are avoided since they have lim- 
ited connection to the underlying physical mechanisms. 
In contrast, mapping closure follows the stretching and 



squeezing of the fluid, and the competition between ram 
pressure, viscosity and advection determines the spectral 
form. 

We study here compressible turbulence without grav- 
ity, self-gravity or thermal conduction. No physical viscos- 
ity is modelled, but numerical viscosity remains present, 
and an artificial viscosity determines the dissipation in re- 
gions of strong convergence. By strong convergence, we 
mean high negative divergence of the velocity field, which 
thus correspond to the shock zones as shown in Fig. 1. 
Periodic boundary conditions were chosen for the finite 
difference ZEUS code simulations, fully described by Mac 
Low ct al. (1998). 

The ZEUS code itself is a time-explicit second-order 
accurate finite difference code (Stone & Norman 1992a, b). 
It is ideal for problems involving supersonic flow and is 
versatile, robust and well-tested. Although higher order 
codes are potentially more accurate, the high speed of the 
algorithms means that large problems can be solved at 
high resolution. This enables us to test for convergence 
of the energy dissipation rate (Mac Low et al. 1998), 
shock distributions (Sect. 2.4) and numerical viscosity 
(Sect. 2.6). Furthermore, the basic hydrocode results have 
been confirmed on solving the same problems with the 
contrasting method of smoothed particle hydrodynamics 
(Mac Low et al. 1998). The constrained transport algo- 
rithm (Evans & Hawley 1988) updated through use of the 
method of characteristics (Hawley & Stone 1995) is used 
to maintain a divergence-free magnetic field to machine 
accuracy and to properly upwind the advection. 

We begin by discussing the method used to count 
shocks from grid-based simulations (Sect. 2.1-2.2). We 
then present the shock jump PDFs and provide analytical 
fits for the hypersonic M = 50 case (Sect. 2.3-2.5). The 
one-dimensional counting procedures are verified through 
a comparison with full three-dimensional integrations of 
the dissipated energy (Sect. 2.6). Supersonic hydrody- 
namic M = 5 (Sect. 3) and magnetohydrodynamic (Alfven 
numbers A = 1 and A = 5) simulations (Sect. 4) are then 
then likewise explored. Note the definition of the Alfven 
Mach number A = v lms /vA, where v\ = B 2 /Aitp where 
u rms is the initial root mean square (RMS) velocity and i>a 
is the Alfven speed. The evolution of the velocity PDFs 
are then presented and modelled (Sect. 5). Finally, we in- 
terpret the results in terms of the dynamical models (Sect. 
6). 

2. Hydrodynamic hypersonic turbulence 

2.1. Model description 

The example we explore in detail is the decay of hyper- 
sonic hydrodynamic turbulence (Fig. 1). The three dimen- 
sional numerical simulation on a D 3 = 256 3 grid with peri- 
odic boundary conditions began with a root mean square 
Mach number of M = 50. The initial density is uniform 
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Fig. 1. Four stages in the early development of the shock field. The distribution of converging regions (i.e. the negative velocity 
divergence of the velocity field) is displayed here for a cross-sectional plane of the hydrodynamic M = 50 simulation. The 
displayed values are the numerical equivalents of the divergence (i.e. the sum of the velocity differences) at each grid location. 
In order of increasing time, the minimum divergence (maximum convergence) and the average divergence within converging 
regions only in the displayed cross-sections are: (-178.4, -46.6), (-441.9, -48.0), (-146.7, -13.74) and (-99.9, -7.3). 



and the initial velocity perturbations were drawn from a 
Gaussian random field, as described by Mac Low et al. 
(1998). The power spectrum of the perturbations is flat 
and limited to the wavenumber range 1 < k < k max with 

k-max 8. 

The simulations employ a box of length 2L (hence L 
is equivalent to 128 zones in this example) and a unit of 
speed, u. The gas is isothermal with a sound speed of 
c s = 0.1m. Hence the sound crossing time is 20L/u. We 
take the time unit as r = L/u. Thus a wave of Mach 50 
would laterally cross the box in a time 0.4r if unimpeded. 

Cross-sections of the shock distribution are shown at 
four times in Fig. 1. We actually display a greyscale im- 
age, biased to the high-convergence regions. An evolution 
is hardly discernable in snapshots showing all converging 



regions equally shaded (as displayed in Fig. 1 of Mac Low 
ct al. 1999). 

2.2. Shock number distribution 

Ideally, we would like to calculate the total surface area for 
each shock strength as a function of time. We introduce 
the shock number distribution function dN/dv, which is 
the number of shock elements per unit shock speed as a 
function of time and shock speed. A shock element is the 
surface area of a shock put into units of the grid cell area. 

To simplify our numerical analysis, we calculate in- 
stead the one dimensional shock jump function. This is the 
number distribution of the total jump in speed across each 
converging region along a specific direction. This is written 
as dN/dvj where Vj is the sum of the (negative) velocity 
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gradients (i.e. E [— 5v x ] across a region being compressed 
in the x-dircction). We employ the jump Mach number in 
the x-direction Mj — Vj/c s rather than Vj since this is the 
parameter relevant to the dynamics. Thus, each bounded 
region of convergence in the x-direction counts as a single 
shock and the total jump in Mj across this region is taken 
as its strength. 

Numerically, over the whole simulation grid (x,y,z), we 
calculate each shock jump through 



(1) 



with the condition that Av x < in the range Xi < x < xj . 
This is then binned as a single shock element. The shock 
number distribution dN/dMj is obviously dimensionless. 

The ID approach neglects both the shock angle and 
full shock strength. The distribution of shock jumps, how- 
ever, is found by adding up an enormous number of con- 
tributions over the whole grid. This method has the ad- 
vantage of being extremely robust, involving no model as- 
sumptions. To be a direct representation of the true shock 
strength function, however, it requires a few assumptions 
to be justified: (1) a one-dimensional shock jump is re- 
lated to the actual shock speed, (2) not too many shocks 
are excluded because their surfaces are aligned with the 
chosen direction, (3) the shock velocities are distributed 
isotropically and (4) unsteepened compressional waves can 
be distinguished from true shocks. 

First, we note that it is an extremely difficult task 
to calculate the actual shock speed for each shock. It 
is, however, unnecessary since the shock speed and one- 
dimensional shock jump are closely related statistically. 
We also take the number of zones at which compressive 
jumps are initiated as the number of shocks (where shocks 
are colliding, the method will be inaccurate). Assumption 
(3) will not hold for the magnctohydrodynamic turbulence 
which has a defined direction. In these cases, the jump dis- 
tributions must be calculated both parallel and transverse 
to the original magnetic field. Assumption (4) will not be 
made: we include all acoustic waves, but we have followed 
the width of the jumps and so can verify whether shocks 
or waves are being counted. This is important since broad 
compressional waves also dissipate energy and become in- 
creasingly significant, of course, as the high-speed shocks 
decay and the flow eventually becomes subsonic. 

Many shock surfaces are distorted, occasionally bow- 
shaped. This does not negate our counting procedure pro- 
vided the curvature is not too strong. Here the relevant 
lengths are the shock radius of curvature and the shock 
width. The latter is determined by our numerical method, 
involving von Neumann & Richtmyer (1950) artificial vis- 
cosity, which here constrains strong shocks to just a few 
zones. As seen from Fig. 1, we can confidently take one- 
dimcnsional cuts across the shock surfaces and equate the 



measured jump to the actual jump in speed of fluid ele- 
ments to a good first approximation. 

In Sect. 2.6, we check our method by demanding con- 
sistency with integrated quantities derived directly from 
the numerical simulations. 

2.3. Hypersonic turbulence 

The random Gaussian field rapidly transforms into a shock 
field in the Mach 50 case (Fig. 1). The shock steepening 
is reflected in the initial increase in the minimum value of 
the divergence (see the caption to Fig. 1). Note that the 
average value of the divergence does not change i.e. the 
total number of converging zones only falls from half to 
about one third, despite the steepening. The explanation 
is that the shocks have time to form in the strongly con- 
verging regions but the compression in most of the flow 
progresses slowly. After the time t = 0.1, the number of 
fast shocks decays and the average convergence decreases. 
The total number of zones with convergence, however, re- 
mains constant throughout. This fact, that the total shock 
surface area is roughly conserved, is verified in the follow- 
ing analysis. 

The one-dimensional distribution of the number of 
shock jumps as a function of time is presented in Fig. 2 for 
the case with RMS Mach number M = 50. This demon- 
strates that the shock jump function both decays and 
steepens. One can contrast this to the decay of incom- 
pressible turbulence where the distribution function, as 
measured by wavenumbcr, maintains the canonical Kol- 
mogorov power-law in the inertial range during the decay 
(Lesieur 1997). Here we remark that a power-law fit is im- 
possible (Fig. 2a) but stress that this result applies only 
to the case at hand: decaying turbulence. 

The jump distribution is very close to being exponen- 
tial in both velocity and time. This remarkably simple 
conclusion is based on the good fits shown in Fig. 2b. Note 
that the pure exponential only applies to the medium and 
strong shock regime. To also account for the low Mach 
number regime, we fitted a further time dependence as 
shown in Fig. 3, yielding 

J^ 10 5 - 72 tcxp(-M^/2.0) ( 2 ) 



in terms of the 1-D jump Mach number. Better fits can 
be obtained with a somewhat more complex time depen- 
dence. We find excellent fits for 

^ = 10 5 - 79 texph/3M^] (3) 

with a ~ 0.88 ± 0.03 and (3 ~ 0.52 ± 0.03. The values 
and errors are derived from parameter fitting of all the 
displayed curves to within a factor of ~ 1.5. We exclude in 
this process the phase where collisional equilibrium would 
not be expected: at early times and low Mach numbers 
tMj < 0.5. Also, we exclude the jump speeds where the 
jump counts are low (shock numbers dN/dMj < 3000). 
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Fig. 2. The jump velocity distribution. Two representations of the number of compressive layers in the x-direction as a function 
of the total x-speed change across the layer at 10 equally-spaced times beginning at time t = 0.5 (solid line) and continuing, 
with a monotonic decay of the high-speed jumps, in time steps of 0.5. They demonstrate the decay and steepening and that (a) 
no power law can represent the behaviour at any stage, and (b) a series of exponentials of the form exp(— t(vj — 2)/2) fit the 
high velocities well. 




Fig. 3. The jump velocity distribution extracted from the M = 50 hydrodynamic simulation as a function of time, as in Fig. 2. 
The fitted function is dN/dM 3 = 10 5 42 t exp(-tM,-/2). 



6 



Smith, Mac Low & Zuev: Supersonic Turbulence 




Fig. 4. The evolution of the jump velocity distribution at two 
resolutions for M = 50 hydrodynamic simulation with kmax = 
2. The shock number for the 64 simulation has been multiplied 
by 4 to adjust for the larger zone sizes. 



Fig. 6. The jump velocity distribution extracted from 
the 128 3 M = 50 hydrodynamic simulation with kmax 
= 2 as a function of time. The fitted function is 
dN/dM, = 10 3 - 93 texp(-tM,-/6). 



2.4- Convergence and dependence on initial conditions 

Resolution studies are essential to confirm numerical re- 
sults. One hopes that the results demonstrate conver- 
gence. This is plausible for supersonic flow in which the 
decay does not depend on the details of the viscosity or 
the details of the shock transitions. This has been con- 
firmed for the analysis of the total energy (Mac Low et al. 
1998). 

We compare available simulations for the hypersonic 
study with 64 3 and 128 3 . We also set the initial wavenum- 
ber range to k max = 2 and can thus examine the depen- 
dence on the chosen initial state. 

The results at the two different resolutions are in 
quite good agreement, especially in the high Mach number 
regime (Fig. 4) . The density of high Mach number shocks 
is quite low and they are well resolved. At low Mach num- 
bers, the lower resolution simulation fails, of course, to 
capture the vast quantities of weak comprcssional waves 
contained in the higher resolution example. It is to be ex- 
pected that shock turbulence configurations get extremely 
complex on small scales, through the interactions which 
produce triple-point and slip-layer structures. To capture 
this structure requires adaptive grid codes. This does not 
mean, however, that the simulations are inaccurate for our 
purposes since energy dissipation is not controlled by the 
weak shocks until very late times, as verified in Fig. 5. 

A similar formula for the shock jump function is found. 
The evolution, however, is three times slower. We show in 
Fig. 6 the model fit 



dN 
dMj 



10 4 - 23 i exp(-M j i/6.0). 



(4) 



This suggests that the rate of decay is proportional to 
the initial mean wavenumber of the wave distribution. The 
decay k ma x — 2 simulation is a factor of three slower than 
in the k max — 8 simulation. The mean wavenumber of 
the two simulations are k m = 1.5 and 4.5, with Hat initial 
energy distributions. Hence, given that Mj = v x /c s and t 
is in units of L/10c s , the shock numbers are approximately 
oc t exp(— k a t Vj) where k Q ~ 1.1 k m . The dependence on 
the initial mean wavenumber is expected since the box 
size should not influence the decay rate if we are indeed, 
as we wish, following the unbounded decay. 

Note that two parameter fitting, as above, in this case 
yields 

^ = 10 4 - 23 tcxp[-/3M^] (5) 
with a ~ 1.02 ± 0.03 and /? - 0.168 ± 0.008. 
2. 5. Shock power distribution 

How is the spectrum of shock jumps related to the decay 
of energy? Here we show that the energy dissipation in 
the fast shocks is directly correlated with their number 
which is decreasing exponentially. Furthermore, the weak- 
est shocks, which merge into an area of 'comprcssional 
waves', are ineffective in the overall dissipation. The result 
is that the moderately-supersonic part of the turbulence 
rapidly becomes and remains responsible for the energy 
dissipation for an extended time. 

The shock power distribution function is here defined 
as the energy dissipated per unit time per unit jump speed 
Vj as a function of jump speed. Here again we employ the 
uni-directional jump Mach number Mj. We actually cal- 
culate the energy dissipated by artificial viscosity acting 
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Fig. 5. The power dissipated as a function of the jump velocity. The power is displayed per unit Mj where Mj = Vj/c s , 
and the data is extracted from the M = 50 hydrodynamic simulation for times t = 1,2,3,4 & 5. The fitted function is 
dE/dMj = 0.016 i Mf 5 exp(-tM,-/2). 



along a specific direction within the shocks as defined by 
convergence along this direction. Hence we anticipate that 
in isotropic turbulence one third of the full loss will be ob- 
tained. The relative contributions of artificial and numer- 
ical viscosity, which also confirms the method employed 
here, are discussed in Sect. 2.5. 

The jump number distribution includes a high propor- 
tion of very weak compressional waves that dissipate lit- 
tle energy. The one-dimensional shock power distribution, 
dE/dMj, shown in Fig. 5, illustrates this. The functional 
fit is guided by the above shock number distributions, 
which we would expect to remain accurate for the high 
Mach number jumps. We display the fit to: 



in Mach number is Mj = M(l- 1/M 2 ). This yields 



dE 
dM~ 



= 0.016 A!?- 5 1 exp(-fc Mj t) 



(6) 



with k = 0.5, which is again remarkably accurate given 
the lack of adjustable parameters. Hence dE/dMj oc 
Mj- 5 dN jdMj. Note that the energy dissipated across an 
isothermal steady shock of Mach number M and pre-shock 
density p is E s = (pc 3 s /2) M 3 (l- 1/M 2 ) and that the jump 



E s = 0.5 pc 3 s M* 



V(l+4/M 2 ) 



(7) 



Therefore, the numerical result suggests a (statistical) in- 
verse correlation between density and shock strength. 

Note that Eq. 6 yields E oc t~ 2 - 5 (on integrating over 
Mj and substituting the variable w = Mjt) Thus, one 
obtains a power-law decay in total energy of the form E 
oc t -1,5 . This behaviour of the total energy decay and the 
energy dissipation rate, is indeed found in the simulations, 
as shown in Fig. 7. Hence the results are fully consistent 
with the simple fits. Note that the energy decay deviates 
from a pure power law: there is an early transition phase 
between times t = 0.1 - 0.5 during which the decay is 
more rapid (Fig. 7). This was not found in simulations of 
moderate Mach number supersonic turbulence (Mac Low 
ct al. 1998) but appears here and also in simulations in 
which the supersonic turbulence is initially driven (Mac 
Low 1999). 



8 



Smith, Mac Low & Zuev: Supersonic Turbulence 




-0.5 0.0 
log(time) 

Fig. 7. The total energy as a function of time (full line) from 
data extracted from the M = 50 hydrodynamic simulation. For 
comparison, the dashed line is a power law with exponent -1.5. 



2. 6. Shock power and artificial viscosity 

In order to check the approximations described in 
Sect. 2.2, we analyse how the energy dissipated due to 
artificial viscosity varies over time and with different pa- 
rameters of the model. For non-driven, or freely decaying 
runs, the following equation is taken as a basic: 



olE 
~dt 



{H v + H numerical) j 



(8) 



where H v and H numerica i are energy dissipation rates due 
to artificial viscosity and numerical viscosity respectively. 

We follow Stone & Norman (1992) Eqs. (32)-(34), 
(134) in computing H u . The first step in computing deriva- 
tives is to realize that artificial viscosity operates on com- 
pression only, so points with positive derivatives in each 
direction are set to zero : 



and compute the artificial viscosity dissipation rate for the 
entire cube at each particular time dump as 



i+l,j,k 



-vl 



if (vl i+ ij,k 
otherwise. 



vli,j,k) > 
(9) 



Then scalar artificial pressures qi in all 3 directions are 
computed, with the uniform mesh (Ax) : 



qi = 1 2 P 







\dx t ) 


_Ax_ 



p(dviY 



(10) 



where (I /Ax) is a dimensionless constant which measures 
the number of zones over which the artificial viscosity will 
spread a shock and was chosen to be 2 in these simulations. 
Then we calculate the artificial viscosity tensor h v : 



h v 



-Vv :Q= (-1) 



dxl 8x2 8x3 



(11) 



H, 



v = J h v dx 3 ~ Y^{h V! ijk)5x 3 . (12) 

ijk 



To understand the convergence properties of the energy 
dissipation rate, we performed a resolution study for res- 
olutions of 64 3 , 128 3 , and 256 3 zones for both hydrody- 
namic and MHD models. Fig. 8 shows that the energy dis- 
sipation rate H v has converged to better than 25% moving 
from the 128 3 to the 256 3 models. 

Now let us examine the fraction of the total kinetic en- 
ergy lost to dissipation in shocks due to artificial viscosity 
R = H v /{dE/dt). We avoid taking this ratio locally in 
time as variations of H v in time around the average lead 
to spurious results. A more robust way to compute this 
ratio is to integrate over the entire curve, thus averaging 
over the variations in H v by computing 



R = 



J Hdt 



Ek(tf) — E k (t ) 



(13) 



Table 1 shows these ratios for a set of models with M = 5 
and kmax = 8. The table indicates that the energy dis- 
sipation ratio R has converged to a few percent in the 
hydrodynamic case at the resolution of 128 3 , and even in 
the MHD case is converged to better than 3% at 256 3 . 

We find that the fraction of energy lost in shocks to 
artificial viscosity doubles when we go from MHD to hy- 
drodynamic models. Mac Low et al. (1998) speculated that 
runs with magnetic fields dissipating much of their energy 
via short- wavelength MHD waves. The factor of two differ- 
ence in the dissipation rate due to artificial viscosity from 
hydrodynamic to MHD runs gives further evidence for 
this dissipation mechanism. This behavior appears well- 
converged, as discussed above. 

3. Supersonic turbulence: M = 5 

Moderate-speed hydrodynamic turbulence has been dis- 
cussed by Mac Low et al. (1998). They also found that 
the kinetic energy decreases with time as a power law, but 
with a shallower exponent. For the M=5 study (Model C) 
they recovered a E oc t -10 law. 

The shock jump distribution for M — 5 is shown in 
Fig. 9. Exponential velocity distributions are again found 
even though the range in Mach numbers over which we 
could expect a specific law is rather narrow (Mj <~ 1 — 3). 
Pure exponential time fits, however, are not accurate. We 
display a suitable fit, of the form 



dN 
dM~ 



= 10 6 - 07 t Q exp(-/3M j r) 



(14) 



with a ~ 0.67±0.05 and /? ~ 0.85±0.1. The corresponding 
shock power distribution is shown in Fig. 10 along with a 
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A 


CO 


CO 
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1 


1 


1 


5 


5 


R 


0.62 


0.68 


0.68 


0.38 


0.36 


0.35 


0.59 


0.49 



Table 1. The fraction of the energy dissipated through artificial viscosity for models of supersonic turbulence. Ratio R is defined 
by Eq. (13). Model labels correspond to those of Table 1 of Mac Low et al. (1998). A is the RMS Alfven Mach number. 




M. = v/c s 



Fig. 9. The jump velocity distribution extracted from the M = 5 hydrodynamic simulation (Mac Low et al. 1998) as a function 
of time. The fitted function is exponential in velocity but more complicated in time (see text) (a — 0.67 and /3 = 0.85 displayed). 
The time sequence shown is t = 1,2,3,5 and 10. 



best-fit family of curves calculated from 

=0.10Mf 6 t Q exp(-/3M^ Q ) (15) 

with a = 0.62 ± 0.04 and (3 = 1.6 ± O.l.Hence, an ex- 
ponential velocity distribution is maintained. The decay, 
however, is slightly slower. Integrating over Mj, yields 
E oc t- 1 - 00 . 

Integrating over Mj , with the limits of integration from 
to oo (to account for all the jumps), and with the sub- 
stitution x — Mt a , we observe that the integral becomes 
time- independent. Thus, integration of Eq. 15 over Mj 



yields dE/dt oc i~ 2 - 23 . This yields dE/dt oc t^ 123 , which 
is quite close to the E oc t^ 1 law found by Mac Low et al. 
(1998). 

Now we want to compare how the two quantities vary 
with time: dE/dt from Eq. 15 and the energy dissipation 
rate due to artificial viscosity H v , derived using the algo- 
rithm discussed in Sect. 2.6. Fig. 11 shows that the two 
methods indeed agree. Thus, the shock power distribution 
method (Sect. 2.6) confirms that the statistical approach 
works with power calculations for each ID converging re- 
gion. 
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Fig. 10. The power dissipated as a function of the jump velocity. The power is displayed per unit Mj where Mj = Vj/c s , and 
the data is extracted from the M = 5 hydrodynamic simulation for times t = 1,2,3,4 & 5. The fitted function, given in the text, 
takes a = 0.62 and (5 = 1.6. 



4. MHD turbulence: M = 5, A = 1 and 5 

An analysis of simulations of MHD turbulence allows us 
to determine which wave modes are involved. The time 
dependence of the kinetic energy of freely decaying MHD 
turbulence has been discussed by Mac Low et al. (1998). 
Remarkably, the kinetic energy also decreases with time 
as a power law, although with only a slightly shallower ex- 
ponent than the equivalent hydrodynamic simulation. We 
consider here the high-field example in which the initial 
RMS Mach number M = 5 and the initial RMS Alfvcn 
number A is unity and the low-field equivalent with A = 

5. Mac Low et al. recovered: E oc t~ - 87 . (at the highest 
resolution of 256 3 ) for the high field case. 

The initial field configuration is simply a uniform field 
aligned with the z-axis. Thus the imposed velocity field 
controls the turbulent energy input, and some energy is 
subsequently transferred into magnetic waves. We impose 
no turbulent diffusion here: magnetic energy may, how- 
ever, still be lost through numerical diffusion or MHD 
wave processes. 

The jump number distribution function for these sim- 
ulation (Fig. 12) possesses exponential velocity distribu- 



tions over a range in Mj. The displayed fit to the high 
field case (Fig. 13) is 

^- = 10 5 - 94 t"exp(-/3M,n. (16) 

where a ~ 0.66 ± 0.03 and (3 ~ 0.62 ± 0.03. For the A = 
5 case, a ~ 0.57 ± 0.03 and ~ 1.01 ± 0.05. 

This analysis indicates that (1) supersonic MHD tur- 
bulence is, mathematically at least, no different from hy- 
drodynamic turbulence in that the shock distribution is 
exponential and (2) the time dependencies are also sim- 
ilar to the hydrodynamic M = 5 case. From Fig. 12, we 
conclude (1) that the distribution of high speed shocks re- 
mains unchanged and isotropic in the low-field case, (2) 
the distribution of transverse waves is somewhat faster to 
decay when a weak field is present, whereas in the strong 
field case (3) the high speed transverse waves survive sig- 
nificantly longer from the outset and (4) the whole spec- 
trum of parallel waves is suppressed by a factor ~ 2. 

The velocity jump distributions plotted here are a com- 
bination of shocks and waves. Due to the high Alfven 
speed, the shocks are predominantly slow shocks in the 
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Fig. 8. Resolution study for 3D models - Energy dissipation 
rate due to artificial viscosity versus time. Top graph: model 
B (64 3 , triangle), model C (128 3 , dashed line), and model 
D (256 3 , solid line). Bottom graph: model N (64 3 , triangle), 
model P (128 3 , dashed line), and model Q (256 3 , solid line). 
For hydrodynamic models we observe that the energy dissipa- 
tion rate H v has converged to better than 35% moving from 
64 3 to the 128 3 models, and to better than 25% moving from 
the 128 3 to the 256 3 models. These values are 37% and 23% for 
the MHD models. Thus, energy dissipation rate due to artificial 
viscosity converges as we go to finer grids. 
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Fig. 11. Top: dE/dt versus time - the graph presents H v for 
model C (grid 128, M=5, hydrodynamic) for 21 time dumps 
(rhombs), and integral of Eq. (8) over dMj (solid line). Both 
are energy dissipation rates due to artificial viscosity. We see a 
close fit of the methods. For comparison, dashed line represents 
total kinetic energy changes dE/dt directly from the numerical 
simulation output. Bottom: the absolute value of the relative 
error between H„ and integrated Eq. (15). 

A = 1 case. These shocks can propagate with wave vec- 
tors in all directions except precisely transverse to the 
field. Their propagation speeds are relatively slow since 
the (initially-uniform) Alfven speed is 5 times the sound 
speed) , and therefore their energy may be transferred into 
faster waves via the turbulence raised in the magnetic 
field. In any case, it appears from Fig. 12 that about two- 
thirds of the compressional wave/shock energy is in trans- 
verse compression modes for the case A — 1. The waves 
in this CclSG, cLS measured here by regions of convergence 
along the axes, are fast magnetosonic waves (close to com- 
pressional Alfven waves). The proportion of each can be 
estimated from the simulations by calculating the jump 
widths (in practice, we here place the extra requirement 
that the one-dimensional jump across each individual zone 
exceeds 0.2 c s , in order to distinguish the shocks from the 
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Fig. 12. The jump velocity distribution extracted from the M 
= 5 low magnetic field (A = 5) and high field (A = 1) MHD 
256 3 simulations (Mac Low et al. 1998) as a function of time. 
The 5 times shown are t = 1,2,3,6 and 10. The two curves 
shown for each time are the velocity jumps transverse to the 
field (dashed) and parallel to the field (solid lines). The fitted 
exponentials are described in the text. 



waves). We then find that at Mj = 1, just over half the 
jumps are slow shocks (with an average resolution of <~ 3.0 
and 3.3 zones transverse and parallel, respectively), while 
for Mj > 4, 76% of the 'jumps' are actually waves (an aver- 
age of ~ 6.0 and 7.7 zones in each converging region). This 
contrasts with the hydrodynamic simulations where, quite 
uniformly, well over 90 % of the jumps are indeed narrow 
shocks, resolved only by the artificial viscosity. These are 
of course estimates which ignore the possibility that many 
flow regions may be quite complex combinations of waves 
and shocks. 

An explanation of why such different types of turbu- 
lence decay in the same functional manner is offered in 
Sect. 6.5. 
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Fig. 13. Exponential fits to the jump velocity distributions for 
the high magnetic field (M = 5, A = 1) MHD 256 3 simulation. 
A mean value for the number of shocks in each direction has 
been taken. The 5 times shown are t = 1,2,3,6 and 10. The 
fitted exponentials are described in the text with a = 0.66 and 
f3 = 0.62. 



5. The probability distribution functions 

A traditional aid to understanding numerically-created 
turbulence is the probability distribution function (PDF) 
of the velocities. Here, we determine the one-dimensional 
Mass Distribution Function by calculating the mass per 
unit Mach number interval of the motion in the x- 
dircction. Note all zones contribute here, whether in the 
shocks or not. 

We recover Gaussian distributions in the velocity, as 
apparent in Fig. 14. This is clearly displayed on a Mass- 
log(M 2 ) plot for the high-speed wings (Fig. 15) where a 
Gaussian would generate a linear relationship. Note that 
each time step produces two lines, one for positive and 
one for negative absolute velocities. At early times, the 
imposed symmetry is still dominant but later on, small 
asymmetries become more apparent. To estimate the time 
dependence we have taken the mean mass fraction of each 
pair, on defining the initial mass density to be unity (i.e. 
a unit mass is initially contained in a box of size L 3 ) and 
found that a fit of the form 



d(mass) 
dM 



= 1.6 t u - lb exp(-0.088 M r ) (17) 



is reasonable (Fig. 16). Interestingly, the decay in time of 
the PDF is not exponential. It is a faster decay law than 
for the shocks. This is inherent to the nature of decaying 
turbulence: in the beginning, closely-following shocks ac- 
celerate the fluid to high speeds. At late times, the fast 
shocks are quite evenly spread out and do not combine to 
produce high acceleration. 
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Fig. 14. The decay of the PDF for the M=50 simulation. The 
distributions are shown for the 8 equal time steps from 0.5 to 
4 



Fig. 16. Fits (dashed lines) to the decay of the PDF for the 
M=50 simulation. The mean distributions (full lines) are shown 
for the 8 equal time steps from 0.5 to 4. The model fits are given 
by Eq. 17, which is tailored so that the total mass in the box 
is conserved. 
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Fig. 15. A log-v 2 display of the decay of the PDF for the M=50 
simulation demonstrates the Gaussian nature of the PDFs. 
The distributions for positive and negative absolute speeds are 
shown for the 8 equal time steps from 0.5 to 4, in descending 
order. 



6. Interpretation of shock number distribution 

6.1. The mapping closure method 

A satisfying understanding of turbulence has been elu- 
sive. One can hope that supersonic turbulence may pos- 
sess some simplifying aspects. We have thus devoted much 
time trying to interpret the above shock distributions of 
decaying turbulence. In this section, we first relate the nu- 
merical simulations to theoretical models which have pre- 



dicted exponential velocity gradient PDFs for other forms 
of turbulence. We then interpret the evolution of the shock 
PDFs with an extension of these models utilising the bal- 
listic (ram pressure) principles behind hypersonic turbu- 
lence. 

We require a model based on local interaction in phys- 
ical space in order to model shock interactions. We adapt 
the heuristic 'mapping closure' model in the form pre- 
sented by Kraichnan (1990), in which analytical approx- 
imations were used to describe the evolution of Burg- 
ers turbulence. First, the initial competition between the 
squeezing and viscous processes is followed. An assumed 
Gaussian reference field is distorted non-lincarly in time 
into a dynamically evolving non-Gaussian field. Velocity 
amplitude and physical space are remapped by choosing 
transformations of particular forms. The PDF of a two- 
point velocity difference changes smoothly from Gaussian 
at very large separations (relating independent points) to 
some function £ at small distances. The mapping functions 
are then determined by matching the evolution equations 
with the dynamical equations. The closure is obtained by 
limiting the form of the distortions to locally determined 
transformations. 

The model for Burgers turbulence provides our inspira- 
tion since, in the hypersonic flow simulations of isothermal 
gas, thermal pressure is only significant within the thin 
shock fronts. Furthermore, individual shock structures are 
predominantly one-dimensional. Care must be exercised, 
however, since regions of vorticity are created behind 
curved shocks which are absent in the one-dimensional 
Burgers turbulence. 
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6. 2. Background formulation 

We assume a known reference field for the initial velocity 
u Q (z) as a function of a reference coordinate system z 
(Gotoh & Kraichnan 1993). The 'surrogate' velocity field 
is u(x, t) and is related to u Q and the reference coordinates 
by vector and tensor mapping functions X and J: 



u, = X l (u ,t), 



dzi/dxj = Jij{\io,£o,t) 



(18) 



(19) 



where reference velocity gradients are £jj j0 = dui^jdzj 
and velocity gradients are &j = dui/dxj. The goal is to 
find the distortions described by X and J for which the 
surrogate velocity field is a valid approximation to the ve- 
locity field as given by the equations of motion. The above 
transformations, however, constrain the allowed forms of 
higher order statistics and, hence, neglect some physics 
which affect the long-term evolution (Gotoh & Kraichnan 
1993). 

Mapping closure also assumes that the velocity PDFs 
P(u ) and P{u) are related to multivariate-Gaussians 
through prescribed forms. The justification is simply the 
statistical mechanics argument, as applied to particle 
speeds in equilibrium thermodynamics. The velocity gra- 
dient PDF is written as Q(£,i) with u and £ taken to be 
statistically independent. 

For convenience, we concentrate on the component of 
the velocity, m, in the direction. We rewrite the decay 
of velocity amplitudes in the simpler form 

u t = X;(u , t) = n(t) u ii0 (zi), (20) 

and the velocity gradients then map through 

dui 



dx. 



so defining Y. The velocity gradient PDF, which contains 
the shock PDF, can now be written 



Si 



^ii,o 



N 
J a 



(22) 



where N(t) normalises the PDF. This is the framework in 
which we can discuss the statistical evolution of velocity 
gradients. 

6.3. Dynamical input 

The momentum equations being solved, with the pres- 
sure gradients neglected, axeDui/Dt = (1/ p){d / dxj)fiaij 
where Dui/Dt = dui/dt + Ujdui/dxj, di.j is the stress 
tensor and \x is the viscosity. Differentiation yields the ve- 
locity gradient equation: 



Dt 



1 d d 

p dxi dxj 



where the usual summation rule applies to j and k (but i 
is a chosen direction). 

To continue, we can derive the evolution of the func- 
tions J, by equating Q(£, t) as derived from these equations 
(yielding a rather complex form of the reduced Liouville 
equation, see Gotoh & Kraichnan 1993) with the Q derived 
from the mapping closure approximation (Eq. 22). Then, 
however, the analysis becomes mathematically dense and 
numerical solutions are probably the best option. 

We here revert to a simple heuristic form of mapping 
closure, taking £ to be any component of £u and the map- 
ping function J = Ja to be determined by requiring the 
probability function Q as derived from substituting for 
.D£ /Dt from the reduced Eq. 23 into the reduced Liouville 
equation, 




Q\=ZQ, 



(24) 



(where [..] C :«,£ denotes the ensemble mean conditional on 
given values of u and £) to be equal to the Q derived from 
the equivalent for mapping closure (see Gotoh & Kraich- 
nan 1993), 



dQ(U) , d (8Y{^ ,t) i 



dt 



+ 



Q&t) =7Q(£,*), (25) 



where 7 = (d/dt) ln(N/J). After some manipulation, this 
yields an equation for the evolution of J of the form (see 
also Eq. 24 of Gotoh & Kraichnan (1993)) 



dJ_ 
~dt 



-r£ J 2 - pk 2 d J 3 + D{ J 3 ) 



(26) 



where fc 2 = ((<9£ /<9z) 2 )/(£ 2 ), angled brackets denoting 
the ensemble mean and D(J 3 ) is a function consisting of 
further non-linear derivative terms of the form J 3 . Also, 
an integral term involving Q is neglected here, which a 
posteriori limits the solutions to high jump Mach numbers 
l£l > r2 Xo {Xo being defined below) 

Note that the left hand side and the first two terms 
on the right are close to the Navier-Stokes form given by 
Kraichnan (1990), with the addition of the function r = r» 
which accounts for compressibility. These terms provide 
the statistics of the field through a non-linear transfor- 
mation of the initial field with known statistics. The fi- 
nal term combines together non-linear derivatives, a step 
which permits further manipulation but with the loss of 
information concerning the constants of integration. 

The evolution of J begins rapidly by the steepening 
of large velocity gradients (first term on the right), until 
balance with the viscous term (second term on the right 
hand side) is achieved. Equating these two terms yields 
the form of the mapping function: J = — r£ /(^/c 2 ). Sub- 
(23) stitution into Eq. 21 yields £ = —r 2 {t)^ 2 /{pkfj. We con- 
vert the initial Gaussian PDF Q (£ ) into <2i(£i) by using 
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Eq. 22 to yield the result that the velocity gradient PDF 
is transformed into an exponential function: 



Qi(6) 



( (O 



N(t) 



■ exp 



16 



2r 2 Xo 



(27) 



where Xo = (£o)/(^d)- This function has been derived 
for the high gradients i.e. the shocks, consistent with the 
numerical simulations. 

The exponential being set up, the gradients are then 
determined by the inviscid terms. Thus, the continued evo- 
lution of J is described by the first two terms in Eq. 26 
which, on substituting the mapping J = £/(r£i) yield 



d_ 

di 



£/6 



(28) 



This has an asymptotic power law solution of the form 
= 1/(1 + t/(ti)), ti being a constant. Hence, at large 
t, the velocity gradient associated with a fluid trajectory is 
inversely proportional to time, as physically plausible for a 
high Mach number expanding flow. Substituting £i = (1 + 
i/ii)£ for £i in Eq. 27, and normalising the PDF through 
N, yields 



Q(0 = fcjexp 



(29) 



for large t, where x — 2 ^o(Co) / (M^dCo) and k is a constant. 
Hence, mapping closure predicts, for zero-pressure hydro- 
dynamic flow, the same fast shock decay law as uncovered 
in the hypersonic simulations. 

The above mapping closure technique provides insight 
into the rapid build up of velocity gradients and transfor- 
mation to an exponential, as well as the evolution of the 
exponential term. The form of the prefactor excludes an 
extension to low speeds. 

6.4. A direct physical model 

The interpretation we now present is an extension of the 
dynamic basis of the mapping closure model. First, we 
neglect viscosity since the long term evolution must be 
derivable from purely inviscid theory. The critical addi- 
tion to the above analysis is a result of the simulations: the 
total number of zones across which the gas is converging 
remains at approximately 30%, independent of time. This 
can also be derived from integrating Eqs. (2) and (4) over 
Mj which yields constant total shock surfaces of 1.05 10 6 
and 0.10 10 6 zone surface elements, respectively (the dif- 
ference being mainly the factor of 8 more zones in the 
former 256 3 calculation). As can be seen from Fig. 3), the 
strong shocks disappear, being replaced by weak shocks. 
We interpret this empirical conservation law as due to 
the fluid being contained predominantly within the layers 
which drive the shocks and the looser-defined layers which 
drive the weaker compressional waves. These layers inter- 
act, conserving the total shock area. This is expected from 



shock theory since the two driving layers merge but the 
leading shock waves are both transmitted or reflected. The 
need for layers to drive the shocks, as opposed to shocks 
to sweep up the layers, is a necessity in an isothermal flow 
where the pressure behind a shock must be associated with 
enhanced density. Nevertheless, a shock will sweep up and 
compress pre-existing density structures ahead of it. 

The decay of a single shock is controlled by the de- 
cay of the momentum of the driving layer. The decay of a 
driving layer is here modelled as due to the time-averaged 
interaction with numerous other layers. These layers can 
be represented by an 'ensemble mean' with the average 
density p. Thus a shock is decelerated by the thrust of 
other shock layers, but its mass is not altered. Mass is 
not accumulated from the oncoming shock layer, but in- 
stead remains associated with the oncoming layer. Then, 
a layer of column density S will experience a deceleration 
of Sdu/dt = -Cdpu 2 where Cd is a drag coefficient of order 
unity and u is the velocity jump (i.e the relative velocity of 
the layer). Integration yields the result ut/L ~ S/(Lp) for 
times exceeding T,/pu where u D is the initial layer speed. 

We impose three physical conditions on the shock dis- 
tribution. 

— At high speeds we take a standard decay law for a 
number of independently-decaying layers. The decay 
rate of fast shocks is proportional to the number of 
shocks present. In this regime, significant numbers of 
new shocks are not generated. 

— Secondly, we shall require that the total number of 
shocks (plus compressional waves, since there is no di- 
viding line in the numerical simulations) remains con- 
stant. The function which satisfies both these condi- 
tions clearly obeys, on integrating over all shocks with 
speed exceeding v\, 



— r 

dt J Vl 



dN 
dv 



dv = —k(vi) 



dN 
dv 



dv 



(30) 



0, to 



dK 
at—e 

dv 



(31) 



with the decay rate function k(v) — > as v — 
conserve the shock number. This has solutions 

dN 

dv 

where a is a constant. 

The third condition we invoke is based on the above 
functional form for individual shock deceleration for 
which vt is a constant. Then, the number of shocks 
above any v = v (t /t) should be conserved i.e. 



dN , 
——dv 

dv 



constant. (32) 

>v (t /t) 

This is similar to the result obtained for the velocity 
gradients in the mapping closure analysis. Integrating 
Eq. 30 with this condition then yields the jump distri- 
bution function 



— = — te 

dv L„ 



vt/L 



(33) 
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where N Q and L Q are constants. 

This implies that we have only two constants with which 
to fit, not a single line, but a whole family of lines! Yet, 
remarkably, this is quite well achieved, as shown in Fig. 3. 
Note that the exponential time-dependence is indeed cor- 
rect at high speeds, but there is a linear time-dependence 
at low speeds, where the weak shocks accumulate. 

6.5. The MHD connection 

We have shown that the decaying MHD turbulence with 
M = 5 possesses similar decay properties to the M — 5 hy- 
drodynamic case despite the different wave phenomena 
involved. Slow shocks, however, possibly dominate the en- 
ergy dissipation in the high-field A=l case. The power 
is dissipated within shock jumps with Mach numbers in 
the range Mj ~ 1 — 2. Alfvcn waves are an important 
ingredient in the exponential tail of the velocity jumps. 
In a uniform medium, Alfvcn waves do not decay even 
when they posess non-linear amplitudes. However, the 
Alfven waves in a turbulent medium will interact non- 
linear ly with other Alfven waves, slow shocks, and den- 
sity structures. Each Alfven wave moves through a mean 
field of other waves. Similarly, each shock layer propagates 
through the 'mean field' of other shocked layers, the basis 
utilised in the above physical model for the hydrodynamic 
case. Hence, the MHD wave interactions could well lead 
to the decay of the velocity jump distribution in the same 
manner as shocks. 

7. Conclusions 

Diffuse gas under various guises is subject to supersonic 
turbulence. Large-scale numerical simulations of 3D MHD 
now allow us to explore many variants. Here, we have stud- 
ied an isothermal gas in which large-scale Gaussian veloc- 
ity perturbations are introduced and freely decay within 
a 'periodic box'. Our main aim here is to analyse the dis- 
tribution of shocks, with the ensuing aim of determining 
the observational signatures. For this purpose, we have 
supplemented the original RMS Mach 5 simulations (Mac 
Low et al. 1998) with hypersonic Mach 50 runs. Indeed, we 
find that the hypersonic case leads to simple mathemati- 
cal descriptions for the shock distribution function, from 
which the M = 5 runs deviate moderately. The Mach 50 
runs obey a steeper energy decay law and (hence) have a 
faster decay of the spectrum of shocks. The velocity PDFs 
remain near-Gaussian but with increasing asymmetry at 
higher speeds and lower masses. 

It should be remarked that the velocity PDFs distri- 
bution and shock distributions decay at different rates. 
The mass fraction at high speeds decreases faster than the 
number of shocks with similar speeds. This indicates that 
the turbulence is indeed decaying from its fully-developed 
state, and the individual shock structures interfere less 



with each other as the flow evolves. That is, the saw teeth 
tend to become more regular with time. 
Further conclusions are as follows. 

— The magnetic field tends to slow down the spectral 
decay as well as the overall energy decay. Fast shocks 
survive longer. 

— Transverse waves of a given strength decay faster than 
waves travelling parallel to the magnetic field. 

— Apart from a small initial period, the energy is not dis- 
sipated by the fast shocks but by the moderate shocks 
with jumps in the Mach number from upstream to 
downstream approximately in the range Mj = 1 — 3, 
even in the M = 50 simulation. This is due to the expo- 
nential fall in fast shock numbers, combined with the 
relatively ineffective dissipation in the weakest shocks. 

Our studies show that for hydrodynamic models the 
ratio of H v to dEk in /dt stays at about 65 percent through 
time, but this ratio varies more with time and has a mean 
value of 30 percent for MHD models. We conclude that 
short wavelength MHD waves are present, and energy loss 
is distributed rather than occuring primarily in thin layers. 

In the hypersonic case, we have found simple laws for 
the evolution and velocity distribution of shock speeds. 
The same basic exponential spectra are also found for 
the very closely related high negative velocity gradients 
of Navier-Stokes and Burgers turbulence, in both simula- 
tions and analytical theory (e.g. Kraichnan 1990, Gotoh & 
Kraichnan 1993). We have extended their heuristic map- 
ping closure model to the present case to reproduce our 
computed spectral forms, and thus demonstrated some of 
the dynamical properties which are inherent to supersonic 
turbulence. 

The relevance of studies of decaying supersonic turbu- 
lence to star-forming clouds was questioned by Mac Low 
et al. (1999). The rapid decay implies that such turbulence 
would be hard to catch in action within long-lived molecu- 
lar clouds. Possible sites, however, within which decaying 
turbulence should prove relevant include the regions down- 
stream of bow shocks, clouds suffering a recent impact and 
disrupted jets. An exponential distribution of shocks such 
as we find will generate very low excitation atomic and 
molecular spectra and inefficient electron scattering, lead- 
ing to steep synchrotron spectra. 

In a following paper, these shock spectra will be em- 
ployed to calculate emission line spectra. On comparison 
with driven turbulence, we may then begin to understand 
the type of turbulence we observe and what may be pro- 
ducing the turbulence. 

Acknowledgements. MDS benefitted greatly from the hospital- 
ity of the Max-Planck-Institut fur Astronomie. We thank E. 
Zweibel for advice and discussions. Computations were per- 
formed at the MPG Rechenzentrum Garching. JMZ thanks 
the American Museum of Natural History for hospitality. Par- 
tial support for this research was provided by the US National 
Science Foundation under grant AST-9800616. 



Smith, Mac Low & Zuev: Supersonic Turbulence 



References 

Evans C, Hawley J.F. 1988, ApJ 332, 659 

Falgarone E., Lis D.C., Phillips T.G. et al. 1994, ApJ 436, 728 

Franco J., Carraminana A., 1999, Interstellar Turbulence, 

CUP, Cambridge 
Galtier S., Politano H., Pouquet A., 1997, Phys. Rev. Lett. 79, 

2807 

Gotoh T., Kraichnan R.H. 1993, Phys. Fluids A, 5, 445 

Hawley J.F., Stone J.M. 1995, Comp. Phys. Comm. 89, 127 

Kraichnan R.H. 1990, Phys. Rev Lett. 65, 575 

Lesieur M., 1997, Turbulence in Fluids, Kluwer (Dordrecht). 

Mac Low M.-M., 1999, ApJ 524, 169 

Mac Low M.-M., Ossenkopf V., 2000, A&A 353, 339 

Mac Low M.-M., Burkert A., Klessen R., Smith M.D. 1998, 

Phys. Rev. Lett. 80, 2754. 
Mac Low M.-M., Smith M.D., Klessen R., Burkert A., 1999, 

Ap&SS 246, 195 
Padoan P., Juvela M., Bally J., Nordlund A., 1998, ApJ 504, 

300 

Porter D.H., Pouquet A., Woodward P.R. 1994, Phys. Fluids 
6, 2133 

Smith M.D., Eisloffel J., Davis C.J. 1998, MNRAS 297, 687 

Stone J.M., Norman M.L. 1992a, ApJS 80, 753 

Stone J.M., Norman M.L. 1992b, ApJS 80, 791 

Stone J.M., Ostriker E.C., Gammie C.F. 1998, ApJ 508, 99 

Vazquez-Semadeni E. 1994, ApJ 423, 681 

Vazquez- Semadeni E., Passot T., Pouquet A. 1996, ApJ 473, 
881 

von Neumann J., Richtmyer R.D. 1950, J. Appl. Phys. 21, 232 



This article was processed by the author using Springer- Verlag 
DTeX A&A style file L-AA version 3. 



